Análisis estadístico

Antes de empezar…

Instala las librerías necesarias (copia y pega en la terminal; no descomentes la línea)…

# install.packages(
#  c("tidyverse", "afex", "emmeans", "GGally", "effectsize", "performance", "see", "qqplotr")
# )

… y carga las librerías más usadas:

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✔ ggplot2 3.3.5     ✔ purrr   0.3.4
## ✔ tibble  3.1.7     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.0
## ✔ readr   2.0.1     ✔ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

El modelado estadístico es difícil, por lo que en este notebook solo cubriremos algunos aspectos básicos. Si en el futuro te enfrentas a experimentos complejos, considera buscar ayuda.

To consult the statistician after an experiment is finished is often merely to ask him to conduct a post mortem examination. He can perhaps say what the experiment died of. (Ronald Fisher)

Modelado estadístico: las dos culturas

En el artículo clásico Statistical modeling: The two cultures encontramos

There are two cultures in the use of statistical modeling to reach conclusions from data. One assumes that the data are generated by a given stochastic data model. The other uses algorithmic models and treats the data mechanism as unknown. (Leo Breiman)

• Explanatory modeling: modelos estadísticos empleados para probar una teoría (causal?). Para ello se emplean variables con un significado científico claro y se evalúa si tienen una relación significativa con la variable de interés. \(\rightarrow\) “Estadística clásica” basada en modelos lineales.
• Predictive modeling: modelos cuyo propósito es predecir observaciones nuevas o futuras. \(\rightarrow\) Aprendizaje automático (machine learning).

En este notebooks nos centraremos en modelos lineales y en los métodos de inferencia clásica. El siguiente notebook cubre modelos basados en aprendizaje automático.

(Otras lecturas interesante sobre las dos culturas es To explain or to predict?.)

Regression is all you need

Regresión simple

El modelo básico sobre el que se construye gran parte de la “Estadística clásica” es el modelo de regresión lineal: \[y = a + b\cdot x + \epsilon\] donde \(\epsilon \sim \mathcal{N}(0, \sigma^2)\).

Es instructivo simular datos que sigan este modelo para entender el significado de la ecuación:

x = seq(-2, 2, 0.1)
expected_behaviour = 2 + 3 * x  # a = 2 y b = 3
epsilon = rnorm(length(x), sd = 1)
y = expected_behaviour + epsilon

df = data.frame('data_x' = x, 'data_y' = y, 
                'expected' = expected_behaviour)
ggplot(df, aes(x = data_x)) + 
  geom_point(aes(x = data_x, y = data_y)) + 
  geom_line(aes(y = expected), col = 2)

# OR...
# plot(x, y)
# lines(x, expected_behaviour, col = 2)

Podemos pensar que expected_behaviour (\(a + b\cdot x\)) es el comportamiento medio o esperado en la población de interés, mientras que epsilon (\(\epsilon\)) representan fluctuaciones aleatorias (bien debidas a errores del proceso de medición o que el proceso estudiado tiene cierta aleatoriedad). Además, la relación entre \(x\) e \(y\) es muy específica. Si aumenta (disminuye) \(x\) aumenta (disminuye) \(y\). Además, un incremento de una unidad en \(x\) siempre produce el mismo incremento en \(y\) (ídem si \(x\) disminuye). Se dice que la relación entre \(x\) e \(y\) es lineal.

La primera pregunta a la que nos enfrentamos es la siguiente. Dados los datos \((x, y)\), ¿podemos estimar expected behaviour?

Ejemplo: linear model (lm)

# 1) crear un modelo lineal
naive_model = lm(data_y ~ data_x, data = df)
# 2) obtener estimaciones de a y b
summary(naive_model)

## 
## Call:
## lm(formula = data_y ~ data_x, data = df)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.79042 -0.61648 -0.07474  0.78419  1.93933 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.0753     0.1439   14.42   <2e-16 ***
## data_x        2.8666     0.1216   23.57   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9215 on 39 degrees of freedom
## Multiple R-squared:  0.9344, Adjusted R-squared:  0.9327 
## F-statistic: 555.5 on 1 and 39 DF,  p-value: < 2.2e-16

# 3) Obtener predicciones del modelo lineal
preds = predict(naive_model, interval = "confidence")
# 4) visualizar el ajuste
data_and_preds = df %>% bind_cols(preds)
ggplot(data_and_preds, aes(x = data_x)) + 
  geom_point(aes(y = data_y)) + 
  geom_line(aes(y = expected, col = "Expected")) + 
  geom_line(aes(y = fit, col = "Predicted")) + 
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.2, fill = "black")

Dado que las estimaciones incorporan el error de estimación, es posible realizar inferencia acerca de la significación de los parámetros.

Ejemplo: inferencia con lm

Un estudiante de biología desea determinar la relación entre temperatura ambiente y frecuencia cardíaca en la rana leopardo, Rana pipiens. Para ello, manipula la temperatura en incrementos de 2ºC que van desde 2ºC a 18ºC, registrando la frecuencia cardíaca (pulsaciones por minuto) en cada intervalo. Los datos están disponibles en “hr.csv”.

# 1) leer los datos
frog_df = read.table("data/hr.csv", header = TRUE)
# 2) Crear un modelo lineal
frog_model = lm(heart_rate ~ temperature, data = frog_df)
# 3) Inferencia
summary(frog_model)

## 
## Call:
## lm(formula = heart_rate ~ temperature, data = frog_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3389 -1.7889 -0.6889  1.7611  5.0111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.1389     1.9170   1.116    0.301    
## temperature   1.7750     0.1703  10.421 1.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 7 degrees of freedom
## Multiple R-squared:  0.9394, Adjusted R-squared:  0.9308 
## F-statistic: 108.6 on 1 and 7 DF,  p-value: 1.629e-05

---

El diseño experimental y los resultados de la inferencia en el ejemplo de las ranas nos invitan a concluirque “el aumento de la temperatura causa un incremento de la frecuencia cardíaca”. Sin embargo, esto no es correcto. Por muy fuerte que parezca la relación entre las variables \(x\) e \(y\), nunca debemos interpretar una variable como la causa de la otra. Una relación significativa entre \(x\) e \(y\) puede ocurrir por varios motivos:

1. \(x\) causa \(y\).
2. \(y\) causa \(x\).
3. Existe un tercer factor (llamado variable de confusión) que, bien directa o indirectamente, causa \(x\) e \(y\).

Por otra parte, respecto a los p-valores…

There is some debate among statisticians and researchers about the appropriateness of P values, and that the term “statistical significance” can be misleading. If you have a small P value, it only means that the effect being tested is unlikely to be explained by chance variation alone, in the context of the current study and the current statistical model underlying the test. If you have a large P value, it only means that the observed effect could plausibly be due to chance alone: it is wrong to conclude that there is no effect (emmeans package authors)

En general, hay consenso en que debe abandonarse la interpretación dicotómica del p-valor (efecto significativo Vs no-significativo, sobre todo teniendo en cuenta que se basan en un threshold arbitrario) y que debe favorecerse los resultados basados en intervalos de confianza y tamaños de los efectos (¡incluso si no son significativos!)

For example, a study on the effects of two different ambient temperatures on paramecium diameter returning an effect size of 20 µm and a p-value of 0.1, if centred on p-value interpretation would conclude ‘no effect’ of temperature, despite the best supported effect size being 20, not 0. An interpretation based on effect size and confidence intervals could, for example, state: ‘Our results suggest that paramecium kept at the lower temperature will be on average 20 µm larger in size, however a difference in size ranging between −4 and 50 µm is also reasonably likely’. (…), the latter approach acknowledges the uncertainty in the estimated effect size while also ensuring that you do not make a false claim either of no effect if p > 0.05, or an overly confident claim. (Lewis G. Halsey, The reign of p-value is over)

Ejemplo: intervalos de confianza en lm

summary(frog_model)

## 
## Call:
## lm(formula = heart_rate ~ temperature, data = frog_df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.3389 -1.7889 -0.6889  1.7611  5.0111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.1389     1.9170   1.116    0.301    
## temperature   1.7750     0.1703  10.421 1.63e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.639 on 7 degrees of freedom
## Multiple R-squared:  0.9394, Adjusted R-squared:  0.9308 
## F-statistic: 108.6 on 1 and 7 DF,  p-value: 1.629e-05

confint(frog_model)

##                 2.5 %   97.5 %
## (Intercept) -2.394015 6.671792
## temperature  1.372241 2.177759

Los resultados sugieren que el ritmo_cardiaco de las ranas se incrementará 1.77 bpms por cada grado Celsius, si bien un aumento de la frecuencia cardiaca en el rango (1.37, 2.17) es igualmente verosímil.

Validación de los modelos

Evidentemente cualquier interpretación está supeditada a que el modelo sea correcto. Debemos ser muy cuidadosos a la hora de verificar que se cumplan las asunciones del modelo de regresión lineal. Podemos usar el acrónimo LINE para recordar las asunciones más importantes del modelo: Linear, Independent, Normal, Equal variances.

Ejemplo: evaluación del modelo naive_model

plot(naive_model, ask = FALSE)

# Comprobar la normalidad con qqplot puede ser difícil. Podemos apoyarnos en 
# performance::check_normality 
library("performance")
# check_normality corre shapiro.test, pero tal y como resalta la documentación
# "this formal test almost always yields significant results for the distribution
# of residuals and visual inspection (e.g. Q-Q plots) are preferable."
is_norm = check_normality(naive_model)

# Para hacer la inspección visual
plot(is_norm)

plot(is_norm, type = "qq")

plot(is_norm, type = "qq", detrend=TRUE)

Ejercicio: predicción de AT abdominal

Després et al. ¹ señalan que el tejido adiposo (AT) se asocia con complicaciones metabólicas consideradas como factores de riesgo para la enfermedad cardiovascular. Es importante, afirman, medir la cantidad de AT intraabdominal como parte de la evaluación del riesgo de enfermedad cardiovascular de un individuo. Després y sus colegas realizaron un estudio para predecir la cantidad de AT abdominal a partir de simples mediciones antropométricas.

Entre las medidas tomadas en cada sujeto, se incluyó la AT abdominal obtenida por TC y la circunferencia de la cintura (datos en Wat.csv’’). Una pregunta de interés es cómo de bien se puede predecir y estimar la AT abdominal a partir del conocimiento de la circunferencia de la cintura. Construye un modelo lineal y plotea sus predicciones. ¿Existe alguna relación entre AT y la circunferencia de la cintura?

at_df = read.table("data/at.csv", header=TRUE, sep = ",") 
names(at_df) = c("subject", "waist", "AT")
at_model = lm(AT ~ waist, data = at_df)
at_with_preds = 
  at_df %>% 
  mutate(fit = predict(at_model))


confint(at_model)

##                   2.5 %     97.5 %
## (Intercept) -259.190053 -172.77292
## waist          2.993689    3.92403

ggplot(at_with_preds, aes(x = waist, y = AT)) + 
  geom_point() + 
  geom_line(aes(y = fit, col = "Predictions")) + 
  xlab("waist circumference (cm)") + 
  ylab("deep abdominal AT (cm2)")

Ejercicio: asunciones en el problema de AT abdominal

¿Se cumplan las asunciones de la regresión lineal en el problema de AT abdominal?

plot(at_model, ask = FALSE)

# No, el ruido es heterocedastico (no es el mismo para todos los valores de x):
# ver plot 3.

Regresión múltiple

Para lidiar con situaciones como la ilustrada en el gráfico de “correlation is not causation” (tiburones Vs helados) necesitamos emplear modelos de regresión múltiple, dado que estos permiten “controlar” las variables de confusión. Crear un modelo de regresión múltiple es análogo al caso unidimensional…

Ejemplo: regresión múltiple

Jansen y Keller[^2] utilizaron la edad y el nivel de educación para predecir la capacidad de dirigir la atención (CDA) en sujetos ancianos. CDA se refiere a los mecanismos neuronales que enfocan la mente en lo que es significativo al tiempo que bloquea las distracciones. El estudio recopiló información sobre 71 mujeres mayores con estado mental normal. En este estudio las puntuaciones más altas se corresponden con un mejor funcionamiento de la atención. Las mediciones en CDA, edad en años y el nivel de educación (años de escolaridad) están recogidos en “cda.csv”. Crea un modelo de regresión múltiple para predecir la puntuación de CDA.

[^2] D. Jansen, M. Keller, Çognitive Function in Community-Dwelling Elderly Women”, Journal of Gerontological Nursing, 29 (2003), 34-43.

library("GGally")

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

df = read.csv("data/cda.csv")
ggpairs(df)

cda_model = lm(CDA ~ AGE + EDLEVEL, df)
# cda_model = lm(CDA ~ ., df)  # equivalente al comando anterior
df$predictions = predict(cda_model)

summary(cda_model)

## 
## Call:
## lm(formula = CDA ~ AGE + EDLEVEL, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.9804 -2.2125 -0.0761  2.2824  9.1230 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.49407    4.44297   1.237 0.220498    
## AGE         -0.18412    0.04851  -3.795 0.000316 ***
## EDLEVEL      0.61078    0.13565   4.503  2.7e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.134 on 68 degrees of freedom
## Multiple R-squared:  0.3706, Adjusted R-squared:  0.3521 
## F-statistic: 20.02 on 2 and 68 DF,  p-value: 1.454e-07

Ejercicio: asunciones del modelo de CDA

plot(cda_model, ask=FALSE)

# todo parece correcto!

---

Hasta ahora solo hemos usado datos continuos, pero nada evita usar datos categóricos como predictores. ¡Ojo! Los coeficientes asociados a datos categóricos no deben interpretarse como una pendiente.

Ejemplo: regresión múltiple con datos categóricos

Construye un modelo de regresión lineal para predecir el peso de una persona a partir de los datos contenidos en “antrop.csv”. Interpreta los coeficientes de la regresión.

antrop = read.csv("data/antrop.csv")
antrop = mutate(antrop, male = factor(male))
antrop_model = lm(weight ~ height + male, data = antrop)
summary(antrop_model)

## 
## Call:
## lm(formula = weight ~ height + male, data = antrop)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5982  -3.5488   0.3738   3.2563  16.3890 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -29.46361    7.31841  -4.026 6.89e-05 ***
## height        0.46974    0.04886   9.615  < 2e-16 ***
## male1         1.23216    0.74242   1.660   0.0978 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.353 on 368 degrees of freedom
## Multiple R-squared:   0.36,  Adjusted R-squared:  0.3565 
## F-statistic: 103.5 on 2 and 368 DF,  p-value: < 2.2e-16

antrop_preds = antrop %>% mutate(fit = predict(antrop_model))
ggplot(antrop_preds, aes(x = height, col=male)) + 
  geom_point(aes(y = weight)) + 
  geom_line(aes(y = fit), lwd = 3)

summary(antrop_model)

## 
## Call:
## lm(formula = weight ~ height + male, data = antrop)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5982  -3.5488   0.3738   3.2563  16.3890 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -29.46361    7.31841  -4.026 6.89e-05 ***
## height        0.46974    0.04886   9.615  < 2e-16 ***
## male1         1.23216    0.74242   1.660   0.0978 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.353 on 368 degrees of freedom
## Multiple R-squared:   0.36,  Adjusted R-squared:  0.3565 
## F-statistic: 103.5 on 2 and 368 DF,  p-value: < 2.2e-16

El modelo se puede escribir como \[weight = -29 + 0.47 * height + 1.23 * is-male\] por lo que 1.23 significa que, de media y para una misma altura, los hombres pesan 1.23 Kg más que las mujeres (hemos ajustado por el efecto de la altura).

Ejercicio: Intervalos de confianza

Usa intervalos de confianza para interpretar los resultados de la regresión.

confint(antrop_model)

##                   2.5 %      97.5 %
## (Intercept) -43.8547524 -15.0724663
## height        0.3736688   0.5658194
## male1        -0.2277644   2.6920896

Aunque es cierto que de media los hombres pesan 1.23 Kg que las mujeres, ¡no podemos descartar que las mujeres pesen más que los hombres! (hasta 0.22 Kg).

Ejercicio: Howell

Los datos contenidos en “howell1.csv” son datos censales parciales del área !Kung San compilados a partir de entrevistas realizadas a finales de la década de 1960. Crea un modelo para predecir el peso de los individuos a partir de la altura y el sexo. Evalúa la bondad del modelo.

howell = read.csv("data/howell1.csv", sep = ";")
howell = mutate(howell, male = factor(male))
howell_model = lm(weight ~ height + male, data = howell)
summary(howell_model)

## 
## Call:
## lm(formula = weight ~ height + male, data = howell)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.4420  -3.6612   0.0892   3.5143  14.2151 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -33.847168   1.093840 -30.943   <2e-16 ***
## height        0.499848   0.007825  63.875   <2e-16 ***
## male1         0.734532   0.432259   1.699   0.0898 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.984 on 541 degrees of freedom
## Multiple R-squared:  0.8858, Adjusted R-squared:  0.8853 
## F-statistic:  2097 on 2 and 541 DF,  p-value: < 2.2e-16

howell_preds = howell %>% mutate(fit = predict(howell_model))
ggplot(howell_preds, aes(x = height, col=male)) + 
  geom_point(aes(y = weight)) + 
  geom_line(aes(y = fit), lwd = 3)

Sin ni siquiera usar plot(howell_model) ya somos capaces de ver que el ajuste es malo… cualquier conclusión basada en un modelo erróneo será errónea (garbage in, garbage out).

Ejemplo: dummy variables y contrastes

El dataset iris (puedes obtenerlo con data(iris)) contiene medidas del sépalo y pétalo de varias especies de iris. Construye un modelo lineal para predecir la longitud del sépalo únicamente en función de la especie. Interpreta los coeficientes de la regresión.

data(iris)
iris_model = lm(Sepal.Length ~ Species, iris)
print(summary(iris_model))

## 
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6880 -0.3285 -0.0060  0.3120  1.3120 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         5.0060     0.0728  68.762  < 2e-16 ***
## Speciesversicolor   0.9300     0.1030   9.033 8.77e-16 ***
## Speciesvirginica    1.5820     0.1030  15.366  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared:  0.6187, Adjusted R-squared:  0.6135 
## F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16

iris_preds = iris %>%  mutate(fit = predict(iris_model))
ggplot(iris_preds, aes(x=Species, fill = Species)) + 
  geom_boxplot(aes(y=Sepal.Length)) + 
  geom_point(aes(y = fit), shape=4, size=3)

Al interpretar los coeficientes de la regresión, observamos que lm ha tomado como referencia la especie setosa. Esto se puede observar usando contrasts.

contrasts(iris$Species)

##            versicolor virginica
## setosa              0         0
## versicolor          1         0
## virginica           0         1

Es decir el modelo es \[sepal = mean-setosa-sepal + 0.93 * is-versicolor + 1.58 * is-virginica.\]

Sin embargo, podríamos reescribir el modelo de otra forma de forma que los coeficientes tengan otro significado. Un ejemplo sencillo sería: \[sepal = mean-versicolor-sepal + \alpha_1 * is-setosa + \alpha_2 * is-virginica.\]

En este caso, simplimente estamos variando la especie de referencia. De hecho, contrast se puede modificar para usar como referencia otro nivel del factor:

Ejemplo: contrastes

contrasts(iris$Species) = cbind("_is_setosa" = c(1, 0, 0), 
                                "_is_virginica" = c(0, 0, 1))
iris_model_2 = lm(Sepal.Length ~ Species, iris)
print(summary(iris_model_2))

## 
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6880 -0.3285 -0.0060  0.3120  1.3120 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            5.9360     0.0728  81.536  < 2e-16 ***
## Species_is_setosa     -0.9300     0.1030  -9.033 8.77e-16 ***
## Species_is_virginica   0.6520     0.1030   6.333 2.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared:  0.6187, Adjusted R-squared:  0.6135 
## F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16

Lo interesante es que podemos ajustar los contrastes de forma que respondan a nuestras preguntas científicas. En general, estos contrastes deben ser ortogonales.

Ejemplo: contrastes ortogonales

Imagínemonos el siguiente universo paralelo. En este universo paralelo solo existe la especie setosa. Una empresa de ingeniería genética te contrata para crear nuevas especies con un sépalo más grande. Desarrollas un método conocido como “Método V”, que tiene dos variantes “V-I” y “V-II”. Los experimentos con estas variantes dan lugar a dos nuevas especies que llamas versicolor (V-I) y virginica (V-II). Te planteas dos preguntas científicas:

1. ¿Es el método V capaz de crear especies con el sépalo más grande?
2. ¿Existe alguna diferencia entre V-I y V-II?

levels(iris$Species)

## [1] "setosa"     "versicolor" "virginica"

contrasts(iris$Species)

##            _is_setosa _is_virginica
## setosa              1             0
## versicolor          0             0
## virginica           0             1

# Usamos versicolor como referencia
contrasts(iris$Species) = cbind(
  "V - setosa" = c(-2, 1, 1),
  "I - II" = c(0, 1, -1)
)
contrasts(iris$Species)

##            V - setosa I - II
## setosa             -2      0
## versicolor          1      1
## virginica           1     -1

v_model = lm(Sepal.Length ~ Species, iris)
summary(v_model)

## 
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6880 -0.3285 -0.0060  0.3120  1.3120 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        5.84333    0.04203 139.020  < 2e-16 ***
## SpeciesV - setosa  0.41867    0.02972  14.086  < 2e-16 ***
## SpeciesI - II     -0.32600    0.05148  -6.333 2.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared:  0.6187, Adjusted R-squared:  0.6135 
## F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16

confint(v_model)

##                        2.5 %     97.5 %
## (Intercept)        5.7602675  5.9263991
## SpeciesV - setosa  0.3599303  0.4774031
## SpeciesI - II     -0.4277344 -0.2242656

El método V parace producir sépalos más grandes. Por otra parte, V-II es mejor que V-I.

---

LLegados a este punto, ¡ya hemos cubierto el 90% de los contenidos habituales de un curso de estadística habitual! Aunque parezca mentira, ya hemos hecho análisis tan complejos como

• Análisis de la varianza (anova): por ejemplo, en el problema del método V,
• Análisis de la covarianza (ancova): con el dataset antrop.csv o howell,
• …

Desde una perspectiva moderna, todos los análisis clásicos (T-test, anova, ancova) pueden considerarse como simples modelos de regresión. Esto demuestra el poder unificador de esta perspectiva. En las siguientes secciones revisaremos sin embargo estos modelos clásicos para afianzar la conexión con los modelos de regresión y profundizar en algunas cosas que nos han quedado en el tintero.

T-tests: comparación de dos grupos

Ejemplo: Comparación de grupos mediante regresión

¿Depende la altura de los !Kung adultos del sexo del inviduo? Emplea los datos en “howell1.csv”.

adult_howell = 
  howell %>% filter(age >= 18) %>% 
  mutate(male = factor(male)) %>% 
  mutate(sex = fct_recode(male, "male" = "1", "female" = "0"))

height_model = lm(height ~ sex, data = adult_howell)

summary(height_model)

## 
## Call:
## lm(formula = height ~ sex, data = adult_howell)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -20.6585  -3.4635   0.2965   3.4715  18.7115 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 149.5135     0.4049  369.25   <2e-16 ***
## sexmale      10.8450     0.5914   18.34   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.537 on 350 degrees of freedom
## Multiple R-squared:   0.49,  Adjusted R-squared:  0.4885 
## F-statistic: 336.3 on 1 and 350 DF,  p-value: < 2.2e-16

Ejemplo: histogramas

Apoya tus conclusiones dibujando el histograma de las alturas por sexo

ggplot(adult_howell, aes(x = height, fill=sex)) + geom_density(alpha=0.2)

Ejemplo: evaluación del modelo

¿Cumple el modelo las asunciones de la regresión lineal?

hist(resid(height_model))

Lo cierto es que, aunque correcto, nuestra aproximación no es la habitual a la hora de comparar la media de dos poblaciones. Si las poblaciones son normales (o si podemos invocar el Teorema Central del Límite), podemos hacer uso del t.test.

Ejemplo: T-test para comparación de medias

t.test(height ~ sex, adult_howell)

## 
##  Welch Two Sample t-test
## 
## data:  height by sex
## t = -18.148, df = 323, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  -12.020596  -9.669319
## sample estimates:
## mean in group female   mean in group male 
##             149.5135             160.3585

Ejemplo: T-test de una sola cola

El fichero “hc.csv” contiene los datos de un experimento para verificar si una hormana de crecimiento funciona o no. ¿Existen diferencias significativas entre el grupo de control y el de tratamiento?

library("effectsize")

## Registered S3 method overwritten by 'parameters':
##   method                         from      
##   format.parameters_distribution datawizard

hc = read.csv("data/hc.csv")
t.test(height_inc_cm ~ group, hc, alternative="less")

## 
##  Welch Two Sample t-test
## 
## data:  height_inc_cm by group
## t = -7.6007, df = 19998, p-value = 1.537e-14
## alternative hypothesis: true difference in means between group control and group treatment is less than 0
## 95 percent confidence interval:
##        -Inf -0.0836217
## sample estimates:
##   mean in group control mean in group treatment 
##                5.989647                6.096364

cohens_d(height_inc_cm ~ group, data = hc, alternative = "less")

## Cohen's d |        95% CI
## -------------------------
## -0.11     | [-Inf, -0.08]
## 
## - Estimated using pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).

# o bien...
height_test = t.test(
  hc %>% filter(group == "control") %>% pull(height_inc_cm),
  hc %>% filter(group == "treatment") %>% pull(height_inc_cm),
  alternative = "less"
) 
cohens_d(height_test) 

## Cohen's d |        95% CI
## -------------------------
## -0.11     | [-Inf, -0.08]
## 
## - Estimated using un-pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).

# ¡o bien podemos dejar que la librería elija el tamaño del efecto adecuado!
effectsize(height_test)

## Cohen's d |        95% CI
## -------------------------
## -0.11     | [-Inf, -0.08]
## 
## - Estimated using un-pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).

Siempre debemos considerar el tamaño del efecto, bien en base a nuestro conocimiento o bien usando los estadísticos adecuados como Cohen’s D. Para Cohen’s D la heurística es:

• \(|d| \approx 0.2\): efecto pequeño (despreciable).
• \(|d| \approx 0.5\): efecto medio.
• \(|d| \approx 0.8\): efecto grande.

Ejercicio: T-test de una sola cola

En los datos de adult_howell, ¿podemos concluir que los hombres son más altos que las mujeres? ¿Cuál es el tamaño del efecto?

t.test(height ~ sex, adult_howell, alternative="less")

## 
##  Welch Two Sample t-test
## 
## data:  height by sex
## t = -18.148, df = 323, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group female and group male is less than 0
## 95 percent confidence interval:
##     -Inf -9.8592
## sample estimates:
## mean in group female   mean in group male 
##             149.5135             160.3585

cohens_d(height ~ sex, data = adult_howell, alternative = "less")

## Cohen's d |        95% CI
## -------------------------
## -1.96     | [-Inf, -1.74]
## 
## - Estimated using pooled SD.
## - One-sided CIs: lower bound fixed at (-Inf).

Ejercicio: T-test apareado

John M. Morton et al. [^3] examinaron la función de la vesı́cula biliar antes y después de la fundoplicatura, una cirugı́a para detener el reflujo. Los autores midieron la funcionalidad de la vesı́cula biliar calculando la fracción de eyección de la vesı́cula biliar (GBEF) antes y después de la fundoplicatura. El objetivo de la fundoplicatura es aumentar la GBEF, que se mide como un porcentaje. ¿Hay evidencia para concluir que la fundoplicatura aumenta el funcionamiento de la GBEF? Datos en “gbef_long.txt” (o “gbef.txt”, para un reto).

# Los datos en gbef.txt no están en el formato adecuado para un análisis estadístico...
# ¡hay que luchar para ordenarlos de forma correcta!

gbef_df = read.table("data/gbef.txt")
gbef_df[3, ] = c("ID", "", 1:(ncol(gbef_df) - 2))
gbef_df = t(gbef_df)

gbef_df = as.data.frame(gbef_df[-(1:2), ]) %>% setNames(gbef_df[1, ])
gbef_df = 
  gbef_df %>% mutate(
  Preop = as.numeric(Preop), 
  Postop = as.numeric(Postop),
  ID = factor(ID)
)

gbef_df_long =
  gbef_df %>%
  pivot_longer(Preop:Postop, names_to="class",  values_to = "gbef")

gbef_df_long = read.table("data/gbef_long.txt", header = TRUE)
gbef_df_long = 
  gbef_df_long %>% 
  arrange(ID)

head(gbef_df_long)

##   ID  class gbef
## 1  1  Preop 22.0
## 2  1 Postop 63.5
## 3  2  Preop 63.3
## 4  2 Postop 91.5
## 5  3  Preop 96.0
## 6  3 Postop 59.0

preop = gbef_df_long %>% filter(class == "Preop") %>% pull(gbef)
postop = gbef_df_long %>% filter(class == "Postop") %>% pull(gbef)

gbef_test = t.test(postop, preop, paired = TRUE, alternative = "greater")
print(gbef_test)

## 
##  Paired t-test
## 
## data:  postop and preop
## t = 1.9159, df = 11, p-value = 0.04086
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  1.131919      Inf
## sample estimates:
## mean of the differences 
##                  18.075

cohens_d(gbef_test)

## Cohen's d |      95% CI
## -----------------------
## 0.55      | [0.03, Inf]
## 
## - One-sided CIs: upper bound fixed at (Inf).

Ejercicio: evaluación del modelo de datos apareados

diff = postop - preop
hist(diff)

ANOVA: comparación de medias para múltiples grupos

ANOVA permite comparar más de dos grupos entre sí (¡como ya hicimos con iris!). Asumiendo los peligros de dar recetas generales, en una primera aproximación podemos seguir los siguientes pasos:

1. Explorar y visualizar los datos.
2. Construir y/o elegir contrastes. ¡Ojo! Si quieres usar sumas de cuadrados de tipo III (recomendado), estos contrastes deben ser ortogonales.
3. Usar el modelo ANOVA y verificar sus asunciones.
4. Calcula contrastes o realiza test post-hoc.

Fíjate que ya hemos cubierto casi todos los pasos en los ejemplos de regresión. Los pasos novedosos son:

• Antes de usar contrastes, usamos un test omnibus (ANOVA). Siguiendo el ejemplo del “método V”, un test omnibus simplemente respondería a “¿Hay diferencias entre la longitud del sépalo entre las distintas especies? El problema de los tests omnibus es que, si existen diferencias, no nos dicen entre qué especies hay diferencias.
• Para descubrir qué especies difieren entre sí podemos emplear: 1) contrastes o 2) tests post-hoc. En general, usaremos contrastes si tenemos hipótesis específicas que deseamos comprobar, mientras que usaremos test post-hoc si no tenemos hipótesis específicas (y queremos descubrir cualquier patrón interesante en los datos).

Ejemplo: el “método V”, otra vez.

Repitamos el análisis realizado para el método V, incorporando además los nuevos pasos.

library("car")  # Anova

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:dplyr':
## 
##     recode

## The following object is masked from 'package:purrr':
## 
##     some

data("iris")

# 1) Visualizar
head(iris)

##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1          5.1         3.5          1.4         0.2  setosa
## 2          4.9         3.0          1.4         0.2  setosa
## 3          4.7         3.2          1.3         0.2  setosa
## 4          4.6         3.1          1.5         0.2  setosa
## 5          5.0         3.6          1.4         0.2  setosa
## 6          5.4         3.9          1.7         0.4  setosa

ggplot(iris, aes(x=Species, y = Sepal.Length, fill=Species)) + geom_boxplot() + 
  coord_flip()

# 2) Especificar contrastes si tenemos alguna hipótesis específica 
# (no te preocupes de los denominadores)
contrasts(iris$Species) =  cbind(
  "V - setosa" = c(-2, 1, 1) / 3,
  "I - II" = c(0, 1, -1) / 2
)

contrasts(iris$Species)

##            V - setosa I - II
## setosa     -0.6666667    0.0
## versicolor  0.3333333    0.5
## virginica   0.3333333   -0.5

v_model = lm(Sepal.Length ~ Species, iris)

# 3) Correr ANOVA: test omnibus
v_model_aov = Anova(v_model, type = 3)

En el caso de modelos ANOVA, existen varios tamaños de efecto. Entre los mas conocidos están eta-squared y omega-squared. Para complicar aún más las cosas, las heurísticas para clasificar el tamaño como pequeño, mediano o largo son distintas que para Cohen’s d. En el caso que nos ocupa, para eta-squared o omega-squared tendríamos:

• Pequeño: \(\approx0.01\).
• Mediano: \(\approx 0.06\).
• grande: \(\approx0.14\).

En general, consulta el siguiente FAQ para consultar las heurísticas.

eta_squared(v_model_aov)  # o effectsize(v_model_aov)

## For one-way between subjects designs, partial eta squared is equivalent to eta squared.
## Returning eta squared.

## # Effect Size for ANOVA
## 
## Parameter | Eta2 |       95% CI
## -------------------------------
## Species   | 0.62 | [0.54, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

omega_squared(v_model_aov) # omega-squared se supone que está menos sesgado que eta_squared

## For one-way between subjects designs, partial omega squared is equivalent to omega squared.
## Returning omega squared.

## # Effect Size for ANOVA
## 
## Parameter | Omega2 |       95% CI
## ---------------------------------
## Species   |   0.61 | [0.53, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

Antes de seguir adelante deberíamos comprobar que se cumplan las asunciones de ANOVA. Dejémoslo para otro ejemplo y sigamos adelante.

En el código anterior, ANOVA nos indica que alguna(s) de las especies tiene(n) un sépalo distinto al de lasotras especies… ¡Sin embargo no nos dice cuáles son estas especies!

Para descubrirlo podemos usar contrastes …

# 4.a) 
summary(v_model)

## 
## Call:
## lm(formula = Sepal.Length ~ Species, data = iris)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.6880 -0.3285 -0.0060  0.3120  1.3120 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        5.84333    0.04203 139.020  < 2e-16 ***
## SpeciesV - setosa  1.25600    0.08916  14.086  < 2e-16 ***
## SpeciesI - II     -0.65200    0.10296  -6.333 2.77e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5148 on 147 degrees of freedom
## Multiple R-squared:  0.6187, Adjusted R-squared:  0.6135 
## F-statistic: 119.3 on 2 and 147 DF,  p-value: < 2.2e-16

confint(v_model)

##                        2.5 %     97.5 %
## (Intercept)        5.7602675  5.9263991
## SpeciesV - setosa  1.0797908  1.4322092
## SpeciesI - II     -0.8554688 -0.4485312

… o usar tests post-hoc. La idea básica de los tests post-hoc es fácil de entender: dado que sé que existe alguna diferencia entre las especies, voy a comparar todas las especies entre sí. El problema es que esto dispara el error tipo I muy rápidamente, por lo que tenemos que ser más conservadores a la hora de aceptar una diferencia como significativa. Esto lo logramos con distintos métodos de ajuste de p-valores. Fíjate que el test omnibus sirve como una primera barrera protectora antes de lanzarnos a hacer comparaciones múltiples.

Entre los métodos de ajuste, podemos distinguir entre 1. Métodos centrados en controlar el family-wise error rate (FWER), cuyo credo es “Si cometo un solo error tipo I todas mis conclusiones se desmoronarán”. 2. Métodos centrados en controlar el false discovery rate (FDR), que se corresponde con el credo (considerablemente más optimista) “vamos a intentar estimar cuántos errores tipo I cometo y a no pasarme de cierto número (pero no pasa nada si hay algún error)”.

Podemos emplear R base para realizar los tests post-hoc…

# Bonferroni es bastante conservador, pero es un ajuste muy conocido
pairwise.t.test(iris$Sepal.Length, iris$Species, p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  iris$Sepal.Length and iris$Species 
## 
##            setosa  versicolor
## versicolor 2.6e-15 -         
## virginica  < 2e-16 8.3e-09   
## 
## P value adjustment method: bonferroni

# Los métodos fdr son "BH" (aka "fdr") and "BY".
pairwise.t.test(iris$Sepal.Length, iris$Species, p.adjust.method = "fdr")

## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  iris$Sepal.Length and iris$Species 
## 
##            setosa  versicolor
## versicolor 1.3e-15 -         
## virginica  < 2e-16 2.8e-09   
## 
## P value adjustment method: fdr

… o bien el paquete emmeans (que tiene ciertas ventajas, como veremos):

library("emmeans")

## 
## Attaching package: 'emmeans'

## The following object is masked from 'package:GGally':
## 
##     pigs

v_model_emms = emmeans(v_model, "Species")
pairs(v_model_emms, adjust="bonferroni", infer=c(TRUE, TRUE))

##  contrast               estimate    SE  df lower.CL upper.CL t.ratio p.value
##  setosa - versicolor      -0.930 0.103 147   -1.179   -0.681  -9.033  <.0001
##  setosa - virginica       -1.582 0.103 147   -1.831   -1.333 -15.366  <.0001
##  versicolor - virginica   -0.652 0.103 147   -0.901   -0.403  -6.333  <.0001
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 3 estimates 
## P value adjustment: bonferroni method for 3 tests

pairs(v_model_emms, adjust="fdr", infer = c(TRUE, TRUE))

##  contrast               estimate    SE  df lower.CL upper.CL t.ratio p.value
##  setosa - versicolor      -0.930 0.103 147   -1.179   -0.681  -9.033  <.0001
##  setosa - virginica       -1.582 0.103 147   -1.831   -1.333 -15.366  <.0001
##  versicolor - virginica   -0.652 0.103 147   -0.901   -0.403  -6.333  <.0001
## 
## Confidence level used: 0.95 
## Conf-level adjustment: bonferroni method for 3 estimates 
## P value adjustment: fdr method for 3 tests

# una de las ventajas de emmeans es que podemos calcular tamaños de efecto 
# para las medias aunque, tal y como se menciona en la documentación: 
# "there is substantial disagreement among practitioners on what is the appropriate 
# sigma to use in computing effect sizes; or, indeed, whether any effect-size 
# measure is appropriate for some situations"
eff_size(v_model_emms, sigma=sigma(v_model), edf = 147)

##  contrast               effect.size    SE  df lower.CL upper.CL
##  setosa - versicolor          -1.81 0.226 147    -2.25   -1.360
##  setosa - virginica           -3.07 0.269 147    -3.60   -2.542
##  versicolor - virginica       -1.27 0.213 147    -1.69   -0.845
## 
## sigma used for effect sizes: 0.5148 
## Confidence level used: 0.95

Ejemplo: Contrastes con emmeans

Una desventaja de codificar los contrastes con contrasts es que hay que prestar atención a los detalles (y el diablo esta en los detalles). Por ejemplo, ¿de dónde sale el 3 del contraste "V - setosa" = c(-2, 1, 1) / 3? Además, cuando los análisis de ANOVA se compliquen las cosas se pondrán realmente feas.

Sirva el ejemplo de contrasts para resaltar que los coeficientes de regresión nos permiten hacer contrastes planeados. Sin embargo, a partir de ahora calcularemos dichos contrastes con emmeans:

emmeans(v_model, "Species") %>% 
  # ¡Fíjate que los denominadores ahora sí son comprensibles!
  contrast(method = list('V-Setosa' = c(-1, 0.5, 0.5), 'I - II' = c(0, 1, -1)), 
           infer = c(TRUE, TRUE)) 

##  contrast estimate     SE  df lower.CL upper.CL t.ratio p.value
##  V-Setosa    1.256 0.0892 147    1.080    1.432  14.086  <.0001
##  I - II     -0.652 0.1030 147   -0.855   -0.449  -6.333  <.0001
## 
## Confidence level used: 0.95

Ejemplo: asunciones de ANOVA

En general, ANOVA asume:

• Las observaciones son independientes dentro de los grupos y entre los grupos.
• Los datos dentro de cada grupo son normales.
• La variabilidad dentro de cada grupo es aproximadamente igual a la
  variabilidad en los otros grupos.

plot(v_model, which = c(1, 2), ask=FALSE)

# o bien
plot(check_normality(v_model), type = "qq", detrend = TRUE)

check_homogeneity(v_model) # oooooohhhhhhh nooooooooooooo :(

## Warning: Variances differ between groups (Bartlett Test, p = 0.000).

Ejemplo: ANCOVA

El dataset anxiety proporciona la puntuación de ansiedad, medida antes y después de aplicar un tratamiento contra la ansiedad, de tres grupos de personas que practican ejercicios físicos en diferentes niveles (grp1: basal, grp2: moderado y grp3: alto). Aunque no tenemos ninguna hipótesis específica, hagamos un análisis de los datos…

anxiety = read.csv("data/anxiety.csv")
head(anxiety)

##   id group pretest posttest
## 1  1  grp1    14.1     14.1
## 2  2  grp1    14.5     14.3
## 3  3  grp1    15.7     14.9
## 4  4  grp1    16.0     15.3
## 5  5  grp1    16.5     15.7
## 6  6  grp1    16.9     16.2

# 1) Vis...
ggplot(anxiety, aes(x = pretest, y = posttest, col = group)) +
  geom_point() + 
  geom_smooth(method="lm")

## `geom_smooth()` using formula 'y ~ x'

# 2) Contrates 
# No tenemos hipótesis :(

# 3) Anova + asunciones
anxiety_lm = lm(posttest ~ pretest + group, anxiety)
plot(anxiety_lm, which = c(1, 2), ask=FALSE)

anxiety_aov = Anova(anxiety_lm, type=3)
print(anxiety_aov)

## Anova Table (Type III tests)
## 
## Response: posttest
##              Sum Sq Df  F value Pr(>F)    
## (Intercept)   0.285  1   1.6808 0.2021    
## pretest     101.289  1 598.3210 <2e-16 ***
## group        74.023  2 218.6293 <2e-16 ***
## Residuals     6.941 41                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(anxiety_aov) 

## # Effect Size for ANOVA (Type III)
## 
## Parameter | Eta2 (partial) |       95% CI
## -----------------------------------------
## pretest   |           0.94 | [0.90, 1.00]
## group     |           0.91 | [0.87, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

# 4) Posthoc tests!
# ...

¡Ojo, queremos comparar diferencias entre las medias ajustadas! La función pairwise.t.test() no usará medias ajustadas, por lo que debemos emplear emmeans:

pairs(
  emmeans(anxiety_lm, "group", adjust = "Tukey"),
  infer = c(TRUE, TRUE)
)

##  contrast    estimate    SE df lower.CL upper.CL t.ratio p.value
##  grp1 - grp2    0.641 0.151 41    0.273     1.01   4.235  0.0004
##  grp1 - grp3    2.985 0.150 41    2.619     3.35  19.862  <.0001
##  grp2 - grp3    2.344 0.151 41    1.976     2.71  15.517  <.0001
## 
## Confidence level used: 0.95 
## Conf-level adjustment: tukey method for comparing a family of 3 estimates 
## P value adjustment: tukey method for comparing a family of 3 estimates

ANOVA factorial

Los diseño de ANOVA factoriales (factorial = más de un factor) permiten el efecto individual y conjunto de uno o más factores. Podemos distinguir varios tipos de análisis factoriales…

• … Diseños independientes,
• Diseños con medidas repetidas,
• Diseños mixtos, …

… que, como veremos, plantean ciertas diferencias en los modelos. En cualquier caso, en los diseños factoriales lo que primero debemos comprender es el concepto de interacción.

ANOVA factorial independiente

Ejemplo: Sin interacción entre los factores principales

Estudiamos el efecto de tres drogas sobre el tiempo de reacción (una de ellas placebo) teniendo en cuenta además el sexo de los pacientes que toman el medicamento. Supongamos que el efecto de las drogas y edad se mide en términos de reducción del tiempo de reacción a algún estímulo y que se obtienen los resultados del fichero “drugs_1.csv”. Visualiza el efecto de las drogas y sexo en los tiempos de reacción y propón un modelo.

drugs_df_1 = 
  read.csv("data/drugs_1.csv") %>%
  mutate(drug = factor(drug), sex = factor(sex))

# interaction.plot(drugs_df$sex, drugs_df$drug, response = drugs_df$response_time)
ggplot(drugs_df_1, aes(x=sex, y=response_time, col=drug)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

drugs_model_1 = lm(response_time ~ sex + drug, data = drugs_df_1)
drugs_df_1$predictions = predict(drugs_model_1)

ggplot(drugs_df_1, aes(x=sex, y=predictions, col=drug)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

Ejemplo: interacciones entre los factores principales

Repite el ejercicio anterior para los datos experimentales “drugs_2.csv”.

drugs_df_2 = 
  read.csv("data/drugs_2.csv") %>%
  mutate(drug = factor(drug), sex = factor(sex))

# interaction.plot(drugs_df$sex, drugs_df$drug, response = drugs_df$response_time)
ggplot(drugs_df_2, aes(x=sex, y=response_time, col=drug)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

drugs_model_2 = lm(response_time ~ sex + drug, data = drugs_df_2)
drugs_df_2$predictions = predict(drugs_model_2)

ggplot(drugs_df_2, aes(x=sex, y=predictions, col=drug)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

Uuuuups, las predicciones son malíiiiiiisimas… El modelo no es capaz de capturar las interacciones presentes en los datos.

Ejemplo: modelado de interacciones

Para modelar interacciones…

drugs_model_2 = lm(response_time ~ sex * drug, data = drugs_df_2) 
# o de forma equivalente
#drugs_model_2 = lm(response_time ~ sex + drug + sex:drug, data = drugs_df_2)
drugs_df_2$predictions = predict(drugs_model_2)

ggplot(drugs_df_2, aes(x=sex, y=predictions, col=drug)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

Ejemplo: completando el análisis para drugs_model_2

Una vez aclarado el concepto de interacción podemos completar el análisis para drugs_df_2.

# 2) Contrastes: aunque hemos dejado contrasts de lado, es recomendable 
# planear los contrastes antes de correr ANOVA. 
# ¡Ojo! ahora creamos listas de contrastes

drugs_contrasts = list(
  "_Drugs-Placebo" = c(1, 1, -2), 
  "_A - B" = c(1, -1, 0)
)

# 3) ANOVA 
drugs_model_2 = lm(response_time ~ sex * drug, data = drugs_df_2)
print(Anova(drugs_model_2, type = 3))

## Anova Table (Type III tests)
## 
## Response: response_time
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 2800.49  1 3985.48 < 2.2e-16 ***
## sex          609.76  1  867.77 < 2.2e-16 ***
## drug        1091.94  2  776.99 < 2.2e-16 ***
## sex:drug     660.92  2  470.29 < 2.2e-16 ***
## Residuals     16.86 24                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

eta_squared(drugs_model_2)

## # Effect Size for ANOVA (Type I)
## 
## Parameter | Eta2 (partial) |       95% CI
## -----------------------------------------
## sex       |           0.71 | [0.53, 1.00]
## drug      |           0.99 | [0.99, 1.00]
## sex:drug  |           0.98 | [0.96, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

# Omitimos los plot de comprobación por sencillez...(comprueba que son correctos)
# plot(drugs_model_2, ask=FALSE)

# 4) contrastes
# Contrastes principales para drogas
emmeans(drugs_model_2, ~ drug) %>% 
  contrast(method = list("drugs" = drugs_contrasts))

## NOTE: Results may be misleading due to involvement in interactions

##  contrast             estimate    SE df t.ratio p.value
##  drugs._Drugs-Placebo   40.132 0.649 24  61.808  <.0001
##  drugs._A - B           -0.997 0.375 24  -2.659  0.0137
## 
## Results are averaged over the levels of: sex

# Contrastes para interacciones 
emmeans(drugs_model_2, ~ drug | sex) %>% 
  contrast(interaction = list("drugs" = drugs_contrasts, "sex" = "consec"), by=NULL)

##  drug_custom    sex_consec    estimate   SE df t.ratio p.value
##  _Drugs-Placebo Male - Female     18.9 1.30 24  14.588  <.0001
##  _A - B         Male - Female     20.2 0.75 24  26.977  <.0001

Lo más interesante del código anterior es reflexionar sobre cada uno de los contrastes y su significado. Los contrastes básicos (sex, Drugs-Placebo, A - B) deberían ser claros pero, ¿qué significan los contrastes para interacciones?

• sex1:drug_Drugs-Placebo: El efecto Drugs-Placebo es diferente en hombres y mujeres? Es decir, ¿el efecto de placebo Vs drogas es comparable en hombres y mujeres?
• sex1:drug_A - B: El efecto drug_A - B es diferente en hombres y mujeres? Es decir, ¿el efecto droga A Vs droga B es comparable en hombres y mujeres?

Ejemplo: visualización de los resultados de un contraste con interacciones

Veamos cómo visualizar que las interacciones entre contrastes y sexo son significativas…

contraste_1 = drugs_df_2 %>% mutate(drug = fct_recode(drug, 'Drug' = 'A', 'Drug' = 'B'))
contraste_2 = drugs_df_2 %>% filter(drug != "Placebo") %>% mutate(drug = fct_drop(drug)) 

ggplot(contraste_1, aes(x=drug, y=response_time, col=sex)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=sex))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

ggplot(contraste_2, aes(x=drug, y=response_time, col=sex)) + 
  stat_summary() + 
  stat_summary(geom='line', aes(group=sex))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

En ambos casos, las pendientes de los grupos son distintas, lo que apoya que hay interacciones en los datos.

Fíjate que, si las interacciones son significativas, la interpretación de los efectos principales no tiene sentido. Por ejemplo, en el contraste 2, la droga B aumenta el tiempo de respuesta para las mujeres, pero lo disminuye para hombres. Por tanto, responder a ¿disminuye la droga B el tiempo de respuesta (para hombres y mujeres)? da una imagen incompleta del problema.

NO INTERPRETES LOS EFECTOS PRINCIPALES SI LA INTERACCIÓN ES SIGNIFICATIVA.

Ejemplo: análisis Posthoc sobre las interacciones

ggplot(drugs_df_2, aes(x=sex, y=response_time, col = drug)) + 
  stat_summary() + 
  stat_summary(geom="line", aes(group = drug))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

drugs_emms = emmeans(drugs_model_2, ~ drug | sex) 
contrast(drugs_emms, interaction = list(drug="pairwise", sex='consec'), by=NULL, 
         adjust = "fdr")

##  drug_pairwise sex_consec    estimate   SE df t.ratio p.value
##  A - B         Male - Female   20.226 0.75 24  26.977  <.0001
##  A - Placebo   Male - Female   19.585 0.75 24  26.122  <.0001
##  B - Placebo   Male - Female   -0.641 0.75 24  -0.855  0.4012
## 
## P value adjustment: fdr method for 3 tests

ANOVA factorial con medidas repetidas

Considera un experimento en el que a cada paciente se le realizan tres mediciones: una al comenzar la prueba, otra medición tras probar el tratamiento A, y otra medición tras probar el tratamiento B…

Los tests F en ANOVA asumen que las medidas son independientes. Cuando se utilizan medidas repetidas como en el ejemplo anterior, se viola esta suposición. Por tanto, cuando tenemos ANOVA con medidas repetidas debemos “cambiar” este supuesto por otro válido. Se asume que la relación entre pares de condiciones experimentales es similar. A este supuesto se le conoce como esfericidad. De forma práctica…

1. La esfericidad se puede evaluar mediante el test de Mauchly.
2. ¿Qué hacer si se viola la asunción de esfericidad? Existen correcciones para intentar aliviar este problema. Las más conocidas se deben a Greenhouse y Geisser (GG) y Huynh y Feldt (HF). La recomendación general es que si la esfericidad es mayor que 0.75 (sobre un máximo de 1) se debe usar HF; en otro caso debe usarse la corrección de GG.

Desafortunadamente, el test car::Anova no incorpora correcciones contra la violación de esfericidad por lo que lo cambiaremos por afex::aov_4.

Ejemplo: ANOVA factorial con medidas repetidas

Existe evidencia de que las actitudes hacia ciertos estímulos se pueden cambiar utilizando imágenes positivas. Como parte de una iniciativa para detener el consumo excesivo de alcohol en los adolescentes, el gobierno financió un estudio para determinar si imágenes negativas podrían usarse para hacer que las actitudes de los adolescentes hacia el alcohol sean más negativas. Cada participante en cada estudio vio un total de nueve anuncios simulados en tres sesiones. En una sesión, cada partipante ve tres anuncios (de cerveza, vino o agua); cada anuncio tiene una actitud asociada (positiva, negativa o neutra). Al cabo de las 3 sesiones, las variables se cruzan completamente produciendo 9 condiciones experimentales. Después de cada anuncio, se pidió a los participantes que calificaran las bebidas en una escala que iba desde -100 (me disgusta mucho) pasando por 0 (neutral) hasta 100 (me gusta mucho). El orden de los anuncios fue aleatorio, al igual que el orden en que las personas participaron en las tres sesiones. Datos en “data/attitude_long.dat” (o “attitude.dat” para un reto).

# Los datos en attitude.dat no son adecuados para anova... hay que arreglarlo!
attitude_df = read.table("data/attitude.dat", header = TRUE)

attitude_df_long = 
  attitude_df %>% 
  pivot_longer(beer_pos:water_neu, names_to = 'name', values_to = 'attitude') %>% 
  separate(name, into=c('drink', 'imagery'), sep = '_') %>% 
  mutate_if(is.character, as.factor)

# 1) plot
attitude_df_long = read.table("data/attitude_long.dat", header=TRUE, 
                              stringsAsFactors = TRUE)
head(attitude_df_long)

##   participant drink imagery attitude
## 1          P1  beer     pos        1
## 2          P1  beer     neg        6
## 3          P1  beer     neu        5
## 4          P1  wine     pos       38
## 5          P1  wine     neg       -5
## 6          P1  wine     neu        4

ggplot(attitude_df_long, aes(x=drink, y=attitude, fill=imagery)) + geom_boxplot()

ggplot(attitude_df_long, aes(x=drink, y=attitude, col=imagery)) + 
  stat_summary() + 
  stat_summary(geom="line", aes(group=imagery))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

library('afex')

## Loading required package: lme4

## Loading required package: Matrix

## 
## Attaching package: 'Matrix'

## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

## ************
## Welcome to afex. For support visit: http://afex.singmann.science/

## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: emmeans() for ANOVA models now uses model = 'multivariate' as default.
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************

## 
## Attaching package: 'afex'

## The following object is masked from 'package:lme4':
## 
##     lmer

# 2) Contrastes: usaremos agua para comparar las bebidas alcohólicas y no-alcohólicas
drinks_contrasts = list(
  "Alcohol-Water" = c(1, -2, 1),
  "Beer-Wine" = c(1, 0, -1)
)
imagery_contrasts = list(
  "Others-Neutral" = c(1, -2, 1),
  "Pos-Neg" = c(-1, 0, 11)
)

# 3) ANOVA y asunciones
attitude_model = aov_4(attitude ~ drink * imagery + (drink * imagery | participant), 
                       data = attitude_df_long, observed=NULL)
# El effect size se corresponde con ges (generalized eta squared)
summary(attitude_model)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                Sum Sq num Df Error SS den Df F value    Pr(>F)    
## (Intercept)   11218.0      1   1920.1     19 111.005 2.255e-09 ***
## drink          2092.3      2   7785.9     38   5.106   0.01086 *  
## imagery       21628.7      2   3352.9     38 122.565 < 2.2e-16 ***
## drink:imagery  2624.4      4   2906.7     76  17.155 4.589e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Mauchly Tests for Sphericity
## 
##               Test statistic p-value
## drink                0.26724 0.00001
## imagery              0.66210 0.02445
## drink:imagery        0.59504 0.43566
## 
## 
## Greenhouse-Geisser and Huynh-Feldt Corrections
##  for Departure from Sphericity
## 
##                GG eps Pr(>F[GG])    
## drink         0.57711    0.02977 *  
## imagery       0.74744  1.757e-13 ***
## drink:imagery 0.79840  1.900e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  HF eps   Pr(>F[HF])
## drink         0.5907442 2.881391e-02
## imagery       0.7968420 3.142804e-14
## drink:imagery 0.9785878 6.809640e-10

# asunciones
plot(check_normality(attitude_model), "qq", detrend = TRUE)

check_sphericity(attitude_model) # ya lo sabíamos :(

## Warning: Sphericity violated for: 
##  - drink (p < .001)
##  - imagery (p = 0.024).

# 4.a) Contrastes ...
emm = emmeans(attitude_model, ~ drink | imagery)
interaction_contrasts = list(
  "drink" = drinks_contrasts, 
  "imagery" = imagery_contrasts 
)

# 4.b) Test post-hoc
contrast(emm, interaction = interaction_contrasts, by = NULL)

##  drink_custom  imagery_custom estimate    SE df t.ratio p.value
##  Alcohol-Water Others-Neutral    -11.4  7.53 19  -1.521  0.1446
##  Beer-Wine     Others-Neutral     15.4  4.77 19   3.242  0.0043
##  Alcohol-Water Pos-Neg           116.8 53.33 19   2.189  0.0413
##  Beer-Wine     Pos-Neg           -63.8 36.02 19  -1.770  0.0928

attitude_emms = emmeans(attitude_model, ~ imagery | drink)
contrast(attitude_emms, interaction = list(imagery = "pairwise", drink = "consec"),
         by=NULL, adjust = "fdr")

##  imagery_pairwise drink_consec estimate   SE df t.ratio p.value
##  neg - neu        water - beer    -6.00 2.31 19  -2.599  0.0265
##  neg - pos        water - beer   -10.00 3.26 19  -3.063  0.0128
##  neu - pos        water - beer    -4.00 3.19 19  -1.255  0.2695
##  neg - neu        wine - water   -12.10 2.33 19  -5.187  0.0003
##  neg - pos        wine - water   -10.75 2.65 19  -4.056  0.0020
##  neu - pos        wine - water     1.35 2.78 19   0.485  0.6334
## 
## P value adjustment: fdr method for 6 tests

# Visualización de apoyo a los contrastes/tests
ggplot(attitude_df_long, aes(x = drink, y = attitude, col = imagery)) + 
  stat_summary() + 
  stat_summary(geom="line", aes(group = imagery))

## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`

Como nota para cerrar ANOVA, una alternativa más flexible a estos modelos son los Linear Mixture Models (LMMs), que pueden manejar datos no-independientes, no-esféricos, etc. Los paquetes principales en R son nlme y lme4.

Ejercicio: Ratas

Como parte de una investigación de agentes tóxicos, se asignaron 48 ratas a 3 venenos (I,II,III) y 4 tratamientos (A,B,C,D). La respuesta fue el tiempo de supervivencia en decenas de horas. Realiza una análisis ANOVA completo de los datos.

Análisis no paramétrico

Hemos visto que la Normalidad de los datos es una asunción importante de los métodos vistos hasta el momento. ¿Qué hacer si los datos violan claramente la normalidad?

Los metodos no paramétricos permiten realizar análisis estadísticos sin asumir ninguna distribución de los datos. A cambio, pierden potencia en los contrastes.

A grandes rasgos tenemos… * Como alternativa a t-test independiente \(\rightarrow\) Mann–Whitney test o Wilcoxon’s rank-sum test \(\rightarrow\) wilcox.test. * Como alternativa a t-test apareado \(\rightarrow\) Wilcoxon signed-rank test \(\rightarrow\) wilcox.test. * Como alternativa a One-way ANOVA \(\rightarrow\) Kruskal–Wallis test \(\rightarrow\) kruskal.test. * Como alternativa a Repeated Measures ANOVA \(\rightarrow\) Friedman’s ANOVA \(\rightarrow\) friedman.test.

Ejemplo: Wilcoxon’s rank-sum

El gasto cardı́aco (litros/minuto) se midió por termodilución en una muestra aleatoria simple de 15 pacientes operados de corazón. Deseamos saber si podemos concluir que la media de la población es diferente de 5.05.

data = c(4.91, 4.10, 6.74, 7.27, 7.42, 7.50, 6.56, 4.64, 5.98, 3.14, 3.23,
5.80,6.17,5.39, 5.77)
w_test = wilcox.test(data, mu = 5.05)
print(w_test)

## 
##  Wilcoxon signed rank exact test
## 
## data:  data
## V = 86, p-value = 0.1514
## alternative hypothesis: true location is not equal to 5.05

w_test %>% rank_biserial()

## r (rank biserial) |        95% CI
## ---------------------------------
## 0.43              | [-0.11, 0.78]
## 
## - Deviation from a difference of 5.05.

Ejercicio: Wilcoxon para datos apareados

En un estudio, se midieron los estreses hemodinámicos en sujetos sometidos a colecistectomı́a laparoscópica. Una variable de interés es el volumen diastólico final ventricular (LVEDV) medido en mililitros. Los datos se encuentra en “lvedv.tsv”. El “baseline” se refiere a una medición tomada 5 minutos después de la inducción de la anestesia, y el término “5 minutes” se refiere a una medición tomada 5 minutos después del “baseline”. ¿Hay cambios en los niveles de LVEDV durante la intervención?

lvedv = read.table("data/lvedv.tsv", header = TRUE)
w_test = wilcox.test(lvedv$baseline, lvedv$X5_minutes, paired = TRUE)

## Warning in wilcox.test.default(lvedv$baseline, lvedv$X5_minutes, paired = TRUE):
## cannot compute exact p-value with ties

print(w_test)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  lvedv$baseline and lvedv$X5_minutes
## V = 11.5, p-value = 0.1139
## alternative hypothesis: true location shift is not equal to 0

w_test %>% rank_biserial()

## r (rank biserial) |        95% CI
## ---------------------------------
## -0.58             | [-0.88, 0.03]

Ejercicio: Wilcoxon de una sola cola

Se ha diseñado un experimento para evaluar los efectos de la inhalación prolongada de óxido de cadmio. Quince animales de laboratorio sirvieron como sujetos experimentales, mientras que 10 animales similares sirvieron como controles. La variable de interés fue el nivel de hemoglobina después del experimento. Los resultados se encuentran en “cadmium.csv”. Deseamos saber si podemos concluir que la inhalación prolongada de óxido de cadmio reduce el nivel de hemoglobina

dat = read.table("data/cadmium.csv", header = TRUE, sep = ",")
w_test = wilcox.test(
  dat %>% filter(class == "exposed") %>% .$hemoglobin,
  dat %>% filter(class == "unexposed") %>% .$hemoglobin,
  alternative = "less"
)

## Warning in wilcox.test.default(dat %>% filter(class == "exposed") %>% .
## $hemoglobin, : cannot compute exact p-value with ties

print(w_test)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  dat %>% filter(class == "exposed") %>% .$hemoglobin and dat %>% filter(class == "unexposed") %>% .$hemoglobin
## W = 25, p-value = 0.003004
## alternative hypothesis: true location shift is less than 0

w_test %>% rank_biserial()

## r (rank biserial) |         95% CI
## ----------------------------------
## -0.67             | [-1.00, -0.39]
## 
## - One-sided CIs: lower bound fixed at (-1).

Ejemplo: Kruskal-Wallis y análisis post-hoc

En PlantGrowth tenemos los resultados de un experimento para comparar los rendimientos (medidos por el peso en seco de las plantas) obtenidos bajo un control y dos tratamientos diferentes para hacer crecer a las plantas.

data("PlantGrowth")
head(PlantGrowth)

##   weight group
## 1   4.17  ctrl
## 2   5.58  ctrl
## 3   5.18  ctrl
## 4   6.11  ctrl
## 5   4.50  ctrl
## 6   4.61  ctrl

# Alternativa: kruskal.test(weight ~ group, data = PlantGrowth)
kw_test = kruskal.test(PlantGrowth$weight, PlantGrowth$group)
kw_test %>% rank_epsilon_squared()

## Epsilon2 (rank) |       95% CI
## ------------------------------
## 0.28            | [0.11, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

pairwise.wilcox.test(PlantGrowth$weight, PlantGrowth$group,
                 p.adjust.method = "BH")

## Warning in wilcox.test.default(xi, xj, paired = paired, ...): cannot compute
## exact p-value with ties

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  PlantGrowth$weight and PlantGrowth$group 
## 
##      ctrl  trt1 
## trt1 0.199 -    
## trt2 0.095 0.027
## 
## P value adjustment method: BH

Ejemplo: Test de Friedman para medidas repetidas

El dataset selfesteem contiene la puntuación de autoestima de 10 personas con obesidad en tres momentos durante una dieta. El objetivo es determinar si su autoestima mejoró con la dieta.

selfesteem = read.csv("data/selfesteem.csv")
head(selfesteem)

##   id time    score
## 1  1   t1 4.005027
## 2  1   t2 5.182286
## 3  1   t3 7.107831
## 4  2   t1 2.558124
## 5  2   t2 6.912915
## 6  2   t3 6.308434

ggplot(selfesteem, aes(x=time, y=score)) + geom_boxplot()

# Alternativa: friedman.test(score ~ time | id, selfesteem)
f_test = friedman.test(selfesteem$score, groups = selfesteem$time, blocks = selfesteem$id)
f_test %>% kendalls_w()

## Kendall's W |       95% CI
## --------------------------
## 0.91        | [0.79, 1.00]
## 
## - One-sided CIs: upper bound fixed at (1).

pairwise.wilcox.test(selfesteem$score, selfesteem$time, paired = TRUE,
                     p.adjust.method = "bonferroni")

## 
##  Pairwise comparisons using Wilcoxon signed rank exact test 
## 
## data:  selfesteem$score and selfesteem$time 
## 
##    t1     t2    
## t2 0.0059 -     
## t3 0.0059 0.0117
## 
## P value adjustment method: bonferroni

Automatización de análisis

Ejemplo: test y plots para varios ficheros

La carpeta “meta” contiene los datos de varios papers de los que se desea hacer un meta-análisis. Como primer paso, deseas verificar que los resultados que obtienes sobre los datos son los mismos que se han reportado en los papers originales. Para ello, leerás los datos (todos tienen dos columnas: response y group), ejecutarás un T-test y generarás un gráfico donde se comparen las respuestas para los niveles de group. Los gráficos resultantes se guardarán en la carpeta “meta/images” y los p-valores de los tests en un solo data.frame p-values.

base_path = "data/meta"
images_path = file.path(base_path, "images")
print(images_path)

## [1] "data/meta/images"

dir.create(images_path, showWarnings = FALSE)
files = list.files(base_path, pattern=".csv")

p_values = data.frame()

for (file in files) {
  df = read.csv(file.path(base_path, file))
  test = t.test(response ~ group, df)
  result = data.frame("filename" = file, "p_value" = test$p.value)
  p_values = rbind(p_values, result)
  
  plot = ggplot(df, aes(x = group, y = response, fill = group)) + geom_boxplot()
  image_name = gsub(".csv", ".png", file)
  ggsave(file.path(images_path, image_name), plot)
}

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Saving 7 x 5 in image

print(p_values)

##        filename      p_value
## 1   paper-1.csv 2.050053e-01
## 2  paper-10.csv 1.170465e-01
## 3  paper-11.csv 3.461138e-07
## 4  paper-12.csv 8.363042e-02
## 5  paper-13.csv 3.838954e-01
## 6  paper-14.csv 2.498351e-05
## 7  paper-15.csv 5.605239e-02
## 8  paper-16.csv 1.036900e-05
## 9  paper-17.csv 6.960640e-01
## 10 paper-18.csv 1.825461e-13
## 11 paper-19.csv 1.300555e-07
## 12  paper-2.csv 5.266617e-01
## 13 paper-20.csv 6.377956e-02
## 14  paper-3.csv 1.030685e-05
## 15  paper-4.csv 2.305727e-04
## 16  paper-5.csv 1.445790e-07
## 17  paper-6.csv 1.603187e-01
## 18  paper-7.csv 1.116642e-09
## 19  paper-8.csv 1.279665e-14
## 20  paper-9.csv 1.254387e-07

write.csv(p_values, "p_values.csv", row.names = FALSE)

---

1. J. Després, D. Homme, M. Pouliot, A. Tremblay, and C. Bouchard, ``Estimation of Deep Abdominal Adipose-Tissue Accumulation from Simple Anthropometric Measurements in Men’’ American Journal of Clinical Nutrition, 54 (1991), 471–477.↩︎